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ABSTRACT 

We present a linear diffusion model for the evolution of the double-peaked 
outburst in the transient source XTE J1118+480. The model treats the two 
outbursts as episodic mass deposition at the outer radius of the disk followed by 
evolution of disk structure according to a diffusion process. We demonstrate that 
light curves with fast-rise, exponential decay profile are a general consequence of 
the diffusion process. Deconvolution of the light curve proves to be feasible and 
gives an input function specifying mass deposition at the outer disk edge as well 
as the total mass of the disk, both as functions of time. The derived evolution 
of total disk mass can be correlated with the observed evolution of the ~ 0.1 Hz 
QPO in the source reported in Wood et al. (2000). 

Subject headings: accretion — black hole physics — binaries: close — stars: 
individual (XTE J1118+480) 

1. Introduction 

X-ray novae (XN), also called Soft X-ray Transients (SXTs), are binary X-ray sources 
that show outburst behavior (Chen et al. 1997). The prevailing theory for their transient 
behavior is based on the disk instability model (DIM) first proposed to explain dwarf novae 
outbursts (see the extensive review by Lasota 2001). Applying this model to XN requires 
significant modifications due to X-ray irradiation from the central disk and the disruption of 
the inner disk, presumably by an advection dominated accretion flow (ADAF), in quiescence. 
The DIM with these modifications is shown to be quite successful in modeling transient 
outbursts similar to the canonical XN A0620— 00 (Dubus et al. 2001). However, it is also 
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clear that some sources undergo outbursts that are very difficult to reconcile with the DIM. 
In this paper we consider one such source, XTE J1118+480, and develop an alternative model 
for the outburst based on modulated mass transfer into the outer disk, rather than a sudden 
increase in viscosity in a pre-existing disk as in the DIM. 

We derive a general prescription for calculating X-ray outburst light curves based on a 
diffusion propagation model in the disk. Wc begin with the same basic diffusion equations 
as Bath & Pringle (1981, hereafter BP81) but wc derive an analytical solution rather than 
solving the equations numerically. We derive a Green's function from the first outburst and 
develop a dcconvolution technique to calculate the mass input rate as a function of time for 
the second outburst. This allows us to derive the time evolution of the mass of the disk 
for the second outburst. Based on the global oscillation model of Titarchuk & Osherovich 
(2000, hereafter TOOO) we expect the QPO to correlate inversely with the disk mass. This 
provides us with at least one observational test for the model. 

It is possible that this model may be apphcable to outbursts observed in other sources 
as well. Brocksopp et al. (2000) point out similarities among five XNe: GRO J0422-I-32, 
XTEJ1118+480, X1354-644, GROJ1719-24, and GS 2023+338. AU five have had at least 
one outburst in which the source remains in the low-hard state through the entire episode. 
In three instances a ~ 0.1 — 0.5 Hz QPO drifting persistently upward in frequency is seen. 
These sources may be good candidates for future tests of this model; however, in this paper 
we limit ourselves to XTE J1118+480 as an example. 

In Section 2 we review the properties of the source, particularly the features of the 2000 
outburst which are difficult to account for in the DIM. In Section 3 we describe the main 
features of the diffusion model. In Section 4 we solve for the Green's funtion for the specific 
case where ijj = 2 and in the general case. In Section 5 we derive the dcconvolution procedure 
to get the accretion rate into the disk and the time evolution of the disk mass during the 
outburst. We apply the full procedure to the data on XTEJ1118+480 in Section 6. In 
Section 7 we show that this model can explain the time evolution of the low-frequency quasi- 
periodic oscillation (QPO) seen during most of the outburst. We conclude in Section 8 with 
a discussion of the successes and limitations of this model. 



2. The X-ray Nova XTE J1118+480 

XTE Jl 118+480 is a black hole (BH) transient that was observed from radio to hard 
X-rays during 2000 January to August. Optical observations after the source returned to 
quiescence confirmed the black hole nature of the compact object by measuring a mass func- 
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tion of /(M) = 6.0± O.36M0 (McClintock et al. 2000; Wagner et al. 2001). This is among 
the highest mass functions measured for any transient. In addition, the combination of the 
small distance to the source (1.8 dz 0.6 kpc; McClintock ct al. 2000) and the high Galactic 
latitude {b = 62°) results in this source having the lowest reddening of any XN. The X-ray 
light curve of the source is peculiar, with an unusual double outburst (see Figure 1). The 
first outburst appears to have a fast rise followed by an exponential decay (FRED morphol- 
ogy), but the second outburst is irregular with multiple maxima. The X-ray spectrum was 
essentially constant as an absorbed power law with photon spectral index of 1.73 (Wood et 
al. 2000), typical of the low/hard state of black holes. A radio counterpart at 5 GHz has 
been detected during outburst, although jets were not directly resolved with MERLIN to a 
hmit of < 65 (d/kpc) AU (Fender et al. 2000). 

An earher paper (Wood et al. 2000) reported detailed evolution of the 0.1 Hz QPO, 
using X-ray data from the USA Experiment and RXTE. The QPO frequency increased 
monotonically from 0.07 to 0.15 Hz over ~ 2 months, failing to correlate with the X-ray 
luminosity which rises and falls over the same interval. In this paper, we explore using 
the disk mass rather than usual quantities (such as luminosity or spectral parameters) to 
correlate with the drifting QPO frequency. The low (~ 0.1 Hz) frequency suggests an origin 
at a large radius in the disk because dynamical timescales scale roughly with the Kepler 
frequency. One theory that could explain such a low-frequency QPO is the global oscillation 
model of TOOO. Their model describes a large-scale coherent oscillation in the outer disk, 
linking QPO issues to those of overall disk structure. The model leads to a BH mass estimate 
for XTEJ1118+480 of ~ IM^ based on the ~ 0.1 Hz QPO frequency, in agreement with 
recent optical data (Wagner et al. 2001) which suggest a BH mass range of 6.0-7.7 Mq. 
(We note that there is a typographical error in equation (14) of TOOO which should have 
x-^^^^ instead of x^^^^^ .) The same QPO was seen in EUV and optical wavelengths (Haswell 
et al. 2000; Hynes et al. 2000). 

The double outburst profile of this source is difficult to account for in the standard DIM. 
If an outburst is triggered by a sudden increase in disk viscosity and leads to the FRED-type 
outburst in 2000 January, how can another outburst, with five times the fluence, occur so 
shortly afterwards? What mechanism turns off the first outburst? And, why does the second 
outburst have such an irregular shape? If the instability is triggered at some maximum value 
of the disk mass, and the amount of mass consumed in the first outburst was replenished in 
~ 30 days, then the mass consumed in the second outburst should be replenished at about 
0.5 years. Thus, the DIM models would have predicted that another outburst should have 
occurred by now. Our alternative model is that the outbursts were fueled by sporadic mass 
transfer from the companion rather than from a large amount of mass stored in the disk 
during quiescence. It is beyond the scope of the current paper to develop the dynamics of 
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the mass transfer which accounts for this behavior. It is instructive to look at the energetics 
of the outburst. The luminosity of the long outburst was about 3 x 10^^ erg/s (for a distance 
of 1.8 kpc) which is only 1.4 x 10~'*LEdd (for an assumed black hole mass of 7Mq), whereas 
the canonical DIM outbursts approach the Eddington limit at their peak. 

Mass transfer modulation models for XN outbursts have previously been considered as 
viable alternatives to the DIM. As the XTE J1118+480 outbursts are not easily reconciled 
with the standard DIM, it seems reasonable to develop a different model for this source. 
For example, Kuulkers (2001) has suggested that the long outburst of XTE J1118+480 is 
a superoutburst of the type seen in the SU UMa class of cataclysmic variables (Warner 
1995), based on striking similarities in the X-ray and optical lightcurves and the appearance 
of superhumps. To account for some of the unusual aspects of XTE J1118-I-480, here we 
investigate an alternative model where the outburst is driven by mass transfer from the 
companion and the diffusive propagation of material in the disk. 



3. Description of the Diffusion Model: Equations and Boundary Conditions 

We now derive the main equation describing diffusive propagation of matter in the disk. 
Using the angular momentum balance [e.g. Shakura & Sunyaev 1973, see also Titarchuk & 
Osherovich 1999 hereafter (T099), eq. 4] 

and the continuity equation 

where M{R, t) and Q{R, t) are the mass accretion in the disk and mass accretion input over 
the disk at the given radius R and at the given time t, and W^p is the r</?— component of a 
viscous stress tensor W. We can get a relation between Wr^ and the surface density E(r, t) 
using equation (5) in T099 with an assumption that the rotational frequency is equal to the 
Keplerian one, cu — cuk' 

Wr^ = -'^VHR^ = (3/2)(GM,)VVEi?-3/2 (3) 

where if is a half-thickness of a disk, G is the gravitational constant, is the mass of the 
central object, rj — pu is a viscosity, p is the density, u — ItVt/S is a kinematic viscosity in 
the disk and It, Vt are the turbulent scale and velocity in the disk respectively. We also use 
the relation E = 2Hp. 
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A relation between M and E can be obtained using equation (1) and (3) with an 
assumption that cu — ujk 

d 

M = 27r X 3R^/^ — [uER^/^] . (4) 
Substitution of M from this equation into equation (2) leads to the diffusion equation 



dt RdR 



Equation (5) is identical to equation (2.3) of BP81. It is worth noting that the derivation 
presented here and BP81 is general, using only the mass and angular momentum conserva- 
tions and an assumption of Kepler rotational velocity in the disk. It does not invoke specific 
assumptions regarding the structure of the disk. The viscosity value, r) = pu implemented 
for calculation of W^p in equation (3) are averaged over the vertical size of the disk. The 
product of 2Hp is in fact an integral of the density, p, over the 2;— coordinate, yielding a 
surface density S. 

In our calculation we will assume that the mass accretion input (presumably from the 
companion) occurs at the outer boundary of the disk at R = Ro, namely 

A{t)S{R - Rq) 

Q{R,t) = ^ . (6) 

Thus, A{t) is the mass input rate as a function of time at the outer edge of the disk. 
Using equations (2) and (6) we derive the differential equation for the disk mass — 
27r J:RdR: 

^ = M(i?o, t) - M{Rin, t) + A{t), (7) 

where mass accretion M{Rin, t) at the inner disk boundary can be calculated using equation 
(4). We assume the reflection boundary at the outer edge of the disk; in other words, the 
net mass flux M{Ro,t) through the outer boundary is zero: 

M(i?o,t) = 0. (8) 

This is compatible with the mass input at the outer edge of the disk from eqn. (6) because 
the mass input is considered to be added just inside the outer radius Rq where the reflection 
boundary condition is imposed. Integration of this equation allows us to get as a function 
of time t: ^ 

Ma{t)^ [ [A{T)-M{R,^,T)]dT. (9) 
^0 
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We replace the function E by the function y{R,t) — R^/'^vTi. Using eq. (4) the mass 
accretion M is then expressed through y as follows 

M = 27r X 3i?^/2|| (10) 
oR 

and the diffusion equation (5) can be written in the form 

dy _ 2>v{R)R''^ d „^dy u{R)A{t)d{R - Rp) 

dt R dR dR 27ri?V2 " ^ ) 

We combine this equation with the boundary condition at the outer boundary 

||-0 at i? = i?o (12) 

which follows from equations (8, 10). We assume that at the inner boundary Rin <^ Rq, 
E = 0, which is equivalent to 

y = at i? = i?in. (13) 
Introduction of a new variable x — R^^^ reduces the main equation (11) to: 

dy_ _ 3u{x)d'^y ^ u{x)A{t)S{x - xq) 

dt ~ " 4x2 dx^ ^ TT ^ ' 

with the boundary conditions 

— — = at X — xo, (15) 
ox 

and 

y = for x^O. (16) 

One can check by direct substitution that the solution of this boundary problem for the 
inhomogeneous equation (14) can presented through the following convolution 

y{x,t)= [ A{T)V{x,t-T)dT (17) 
Jo 

where V satisfies the homogeneous equation 

^ = ^-^=^^ (18) 



with the initial condition 



^^^^ ^ K^)5(x-Xo) ^^9^ 

TT 
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and boundary conditions (15-16). Because our ultimate goal is to find the X-ray luminosity, 
which is presumably related to the energy release of the accreting matter at the inner disk 
edge, wc will calculate M at R = -Rin (or at a; — 0). We assume that the mass accretion 
rate at the inner edge is converted with efficiency £eff into the X-ray luminosity, i.e. Lx{t) — 
£effM(0,t). Thus 

L,{t) = £efrM(0, t) oc M(0, t) = = 3n f A(r) V ~ ^r. (20) 



4. The analytical solution of the main problem: General and particular cases 

The solution V{x, t) of equation (18) with the initial condition (19) at i = and bound- 
ary conditions (15-16) can been presented using separation of variables as a series 

T/r^ +\ _ ^-Xlt ^n{x)Xn{xQ)p{xQ) u{xq) 



n=l 



where Xn{x) and A„ are eigenfunctions and eigenvalues which can be found from the inho- 
mogeneous ordinary differential equation: 

X: + Xlp{x)Xr^ = (22) 

combined with the boundary conditions 

Xn = for x^ 0, (23) 

^ = at x = Xo. (24) 

dx 

I I is the norm of the eigenfunction, which is calculated through the integral (for example, 
see the derivation of this formula in Titarchuk, Mastichiadis, & Kylafis 1997) 

||X„||2= / p{x)Xl{x)dx, (25) 
Jo 

where p{x) — 4x^/3i/(x) is the weight function. 

Here, we consider a general class of problems where i/(a;) = i>oX^. In the following 
sections, we solve the problem for one specific value of ■0, then present the general solution. 
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4.1. Case with ip — 2 



When the turbulent scale It is proportional to the characteristic scale R and the turbulent 
velocity Vt is independent of radius, ip = 2 (i.e. u = uqR). In this case, 



Xn{x) — sin 



L(3z/o)V2 

and the eigenvalues A„ are found from the equation 

2 



cos 



where the solutions of this equation are 



An — 



(3i/o)'/'7r(2n-l) 



(26) 



(27) 



(28) 



2 2xo 

For this particular case the solution is expressed by harmonics (eqs. [21] and [26]) 
V{x, t) = V exp[-7r2(2n - 1)^4^0] sin[7r(2n - l)x/2xo] sin[(2n - l)7r/2] (29) 

n=l 

where to = 4a;o/3i/o = 4i?Q/3z/(i?o)- is the viscous timescale and determines both the rise 
and fall time of the response function. 

We define the Green's function K(t) as the flux through the inner edge of the disk at 
time t from delta-function injection at the outer edge at time t = 0. Thus, it is proportional 
to dV{0,t)/dx: 



K{t) oc 



(9^(0, t) 
dx 



i/oXlexp[-7r^(2n - lft/AtQ]{2n - 1) sin[7r(2n - l)/2]. (30) 



n=l 



By definition, K{t) should be normalized to unity, which can be done analytically in the case 
ip — 2. The integral of the series over t in this equation is X]^i(— l)"^~^(2n — 1)~^ = to/n. 

The asymptotic form of the series (29) for t <^ 4io/7r^ can be presented as the integral 



V{x,t) 



2xlvo 



exp[— A Xoi/^o] sin(Aa;) sin(Axo)ciA, 



(31) 



where A = 7r(2n — l)/2a;o and consequently d\ = tt/xq. Transformation of the product of 
sin(Aa;) sin(Aa;o) into the difference [cos A(a; — Xq) — cos A(a; + Xq)\/2 followed by integration 
over A leads us to the formula 

2 1/2 

V{x, t) c ^^-^-j-^ {exp[-(x - xof/ixl{t/to)] - exp[-(x + xof /ixl{t/to)]] . (32) 
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After differentiation we have 

K{t) oc ^^^M = ^3 X 2[^(^y^^)]3/2 exp(-to/4t) for t « 4to/7r^ (33) 

At large times, we can use the series (30) to get 

K{t) oc ^XpH = ^^x exp[-7r't/4to] for t > 4to/7r\ (34) 

From a combination of the two asymptotic forms (eq. [31,34]) wc construct the simple 
analytical approximation of the Green's function. This combination is a product of the 
combined exponentials with a linear superposition of the factors before the exponents with 
weights 1 and 2 respectively: 

K{t) ~ [1 + (to/vrt)^/'] exp[-(to/4t + 7:H/4to)], (35) 

where C is a normalization constant to make the integral of K{t) equal to unity. We have cho- 
sen the weights of the two terms in this approximation to make it more accurately reproduce 
the series result. 



4.2. General Case 

In this section we investigate the asymptotic behavior of the Green's function for t <^to 
and t ^ to in the case of arbitrary kinematic viscosity as a function of R {or x — y/R). We 
have already demonstrated for the ■0 = 2 case that the asymptotic for t <^to is determined 

by the contribution of terms of the series with n 2> 1 (see eq. [30]). The calculation of the 
series is reduced to an integral over A (eq. [31]) where the contribution of terms with n of 
order few is negligible. On the other hand, for t ^ to the term with n = 1 is dominant and 
the contribution of terms with n > 1 is negligible. It is worth noting that the integrated 
function in equation (31) is a product of an exponent and harmonics. The similar asymptotic 
form can be obtained in a general case of the kinematic viscosity. 

One can check by direct substitution that the eigenfunctions X„(x) defined as nontrivial 
solutions of equations (22-24) have the harmonic asymptotic form for large n 

X„{x) = p~^^'^{x) cos[A„/x(x) + if] (36) 

with /x'(x) = dji/dx = ^/p{x) when 

A„ » max{|(pV2)Y^i/2|^ i^y^i^ [\{p-'^r Ml'^'}- (37) 
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They form an orthogonal set of functions with the squared norm 

||X„||^= / p{x)Xl{x)dx = / jj!{x)cos^\n{x) + ^]dx:^ n{xo)/2. (38) 
Jo Jq 

With an assumption that the kinematic viscosity z/ is a power law function of the radial 
coordinate x, namely v = ^qx"^, one can find that Xn is expressed through the Bessel 
functions 

X„ = C„x^/Vi/(4_^)[/i(a;)A„] (39) 

with ii(x) = 2/{A — ip) X (4/3z/o)"^^^3;(^~'^)/^. Among the two Bessel functions Ji/(4_^) and 
J_i/(4_^), only the former one allows Xn{x) to satisfy the inner boundary condition X„ — > 
for x ^ 0. 

We now investigate the behavior of the solution X„ in the limits of large and small 
arguments n{x)\n- Using the asymptotic form of the Bessel function for large arguments 
(Abramowitz & Stegun 1968) 



'K^i{x)\r. 



1/2 



cos[A„/i(x) - 7r(6 - V)/4(4 - V)]. (40) 



Comparing this formula with equation (36) wc sec that Cn = [7rA„/(4 - V)]^/^- For small 
values of the argument, the eigenfunction X^ can be presented as series expansion around 0, 

^'\^'^\uix)\ /2l'/(^-'^) V [-(An/^(3^))V4]^ 



XJx) 



k\r[{5 - V')/(4 -i^) + kY 



.(4-V'). 

In fact, it follows from this equation that 

Xn{x) oc X (42) 

when Xnfi{x) -C 1. 

The eigenvalues A„ are determined by the equation 

^ = for x^xo. (43) 

ax 

For small n = 1,2 the eigenvalues are best determined using the scries (41), whereas for 
large n the asymptotic form (40) is appropriate. Fortunately, even Ai calculated using the 
two different methods produces results which differ from each other by less than one percent. 
For example, the first three terms in series (41) provide an accuracy of Ai better than 0.1%. 
The use of these three terms leads us to a quadratic equation ior u — {A — ip)z'^ /A where 
z = X„ii{xo): 

-2{5-ip)u + 2{5-ip) ^0. (44) 
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The first root of this equation 

= 5 -^-[(5-^^)2- 2(5 (45) 

allows us to calculate the first eigenvalue. For example, for = 1, 2^1 = 1.25, whereas using 
the second method gives zi — 1.18 (see below). 

For calculation of the eigenvalues using the eigenfunction form (40) we have the tran- 
scendental equation: 

tan[z - 7r(6 - V')/4(4 - ^)] = -{2 - ^)/2z{A - ^p). (46) 

The solution of this equation can be written in an analytical form: 

z^[{n- l)7r + 7r(6 - V')/4(4 - ^P) - (47) 

where 

(48) 



2(4-'0)[(n- l)7r + 7r(6 --0)74(4 -V')]' 
We apply equations (47-48) for calculations of the eigenvalues for any n because even for the 
first eigenvalue we have almost the same results using these two methods (compare with Eq. 
45). 

The Green's function K{t) can be written using (41) as follows: 

K{t) oc 1^(0, t)^D'Y, e-^"*A^X„(xo), (49) 

n=l 

where 7 = (6 — ■0)/2(4 — -0) and 

- 'Slh (3^)"'"*' - - ^)/(4 - n- (50) 

Using the asymptotic form of aX, x — Xq (see eq. 40) we can rewrite equation (49) 
as follows: 

dV °° 

K{t) oc -^(0,t) ^Dx ^A>xp(-A^i)cos(6A„ -7r7/2), (51) 

n=l 

where h — t^^^ and D — D'(3uo/ Ay/^x^'^^^'*^ . As in the case with ijj — 2, the viscous time is 
to which in the general case is 

*° = "'(^"^ = 31;^ (4^^- (^2) 
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Thus for t > = (zf/to)-^ or t > to/zf 
dV 

— {Q,t)^Dx {zl/toy/^ cos[2(l - 7)/7r(7 - 3)] ew{-zft/to). (53) 

For t <S io/-2^i one should take into account all terms in the series (eq. [51]). As we have 
already demonstrated in §3, the summation of such a series can be reduced to an integration 
over A: 

dV f^^"^ r°° 

—(0, i) ~ L» X ^ / A^ exp(-A2i) cos(6A - n^/2)dX, (54) 
ax 71 Jo 

1 /2 

where dA ~ tt/^q . 

Now we proceed with calculation of the integral from equation (54), 

/•oo 

7 = / A^exp(-A2t)cos(6A-7r7/2)ciA- cos(7r7/2)/c + sin(7r7/2)/^, (55) 

where 



Ic= / A^exp(-A"t)cos(6A)dA (56) 
Jo 

and 

poo 

Is= X'^exp{-XH)sm{bX)dX. (57) 
Jo 

Expansion of cosine and sine in series and integration over A leads us to (see also Prudnikov, 
Bruchkov & Marichev 1981) 

Ic = ^i-(^+^^/'r[(7 + 1)/2]M[(7 + l)/2, 1/2, -by At] (58) 

and 

7, = ^rM/2r[(7 + 2)/2]M[(7 + 2)/2, 3/2, -b'/it]. (59) 

The asymptotic form of the degenerate hypergeometric function M(a, b, z) = e^z"'~''r{b)/r[a) 
for large arguments z (Abramowitz & Stegun 1968) gives us the asymptotic form for integrals 
Ic and Ig-. 

{/e, /.} ^ (7rVV2)(io/4)-(^+i)/2(to/4t)^+^/' eM-to/^t) for t « to/4. (60) 
Finally we have 

— (0,i) ~ Dto^/^ X7r-^/^(V2t)^+^/^exp(-to/4t)sin(7r7/2 + 7r/4) when t<^to/A. (61) 
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We construct the analytical form of the Green's function K{t) using the asymptotic forms 
cq. (53) for t ^ to and eq. (61) for t <S to, similar to what was done for ^0 = 2 in Section 



ir(t) = C-H2.57r-i/2(to/2t)^+^/'sin(7r7/2+7r/4) + (^2)7/2co3p(l-7)/7r(7-3)]}exp(-to/4t-Zi2tAo), 

(62) 



where C is a normalization constant to make the integral of K{t) equal to unity. The factor 
2.5 in the first term was empirically determined to make this combination of asymptotic 
forms a better match to the series formulation. This form is primarily a computational 
convenience. Note that this form breaks down for ■(/' > 3 so we restrict ip to be less than 
3. A graphical comparison of the series and analytical forms for several values of ip can be 
found in Figure 2. 

An important consequence of the previous derivation is that the FRED outburst shape 
is a natural consequence of the diffusion model. The exponential decay is due to diffusion 
in a bounded medium, while the fast rise is due to the fact that the probe (X-ray emission 
from the inner edge of the disk) is located at a distance from the source (accreting material 
added to the outer edge of the disk) . Thus we have derived the outburst shape for a large 
class of viscosity structures when the viscosity is a function of only the radius. 



In this section, we demonstrate a deconvolution procedure to derive A{t), the mass 
input rate at the outer edge of the disk, and then use A{t) and Lx{t) to derive the mass of 
the disk as a function of time. Now we proceed with general case of the input mass supply 
in the disk Q{R,t) which is proportional to A{t) (see eq. 6). Furthermore, we assume that 
the mass accretion rate at the inner disk edge is converted with efficiency £eff into the X-ray 
luminosity, i.e. Lx{t) — eesM{t, Rm)- Thus using the light curve of Lx{t) we can restore 
information regarding the input function A{t) by the deconvolution of the following equation 
(see eq. 20): 



To do this we also assume that the variability time scale for A(t) is much longer than the 
viscous time scale tg. In other words the support T of the function A{t) is much larger than 
that for K{t). Then using the steepest descent method we get from equation (63) that 



4.1: 



5. Solving The Inverse Problem 




(63) 




for t < T, 



(64) 
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or ^ 

L^{t) = eM{0, t) ~ 37r£eff>l(w) ^ K{t - T)dT for f > T, (65) 

where Tmax is the time point for the maximum of K{t — r) in the integration interval. The 
Green's function K{t) has a sharp maximum at (which depends on tjj and to) and thus 
Tmax = for < t < t*, T^ax = t - ioT < t < T , and Tmax = T for T < t < T + 1*. With 
the steepest descent approximation of the convolution (64-65) one can resolve that equation 
with respect to the input function A{t): 

A{t - tos) oc L^{t)/ [ K{t')dt' for t < T, (66) 
Jo 

where tos = min(t,t*), and 

A{t-tos)ocL^{t)/ [ K{t')dt' for t > T, (67) 

Jt-T 

where t^s — max[(i — T),t^. Notice that ^off = ^ — ^max- 

To calculate the mass of the disk as a function of time, we make the assumption that 
the only terms that matter are mass input at the outer boundary of the disk (proportional 
to A{t)) and the mass leaving the disk into the black hole (proportional to L^). Hence the 
instantaneous change of the disk mass Md{t) is determined by the difference between the 
input function A{t) and the output function. If we assume that the output mass accretion 
rate M(0,t) is proportional to the observed luminosity Lx{t) in accordance with equation 
(20) and the input function A[t) is found as a result of deconvolution with the operator K"-*- 
defined by eqs. (66-67), we can calculate the disk mass as a function of time as follows: 

Md(t)oc / {^-\L,{t)]-L,{t')}dt'. (68) 

It is seen from eqs. (66-67) that K"-*^ is a composite operator which shifts and differentially 
expands the input function Lx{t). 

6. Application to XTEJ1118-I-480 

We now apply this formalism to the 2000 January-August outburst of XTE J1118+480. 
Figure 1 shows a merged 2-12 keV hght curve from the ASM on RXTE and the USA 
Experiment on ARGOS, which together cover the full double-peaked outburst. Only the 
ASM viewed the initial (FRED) outburst because the notification trigger did not come 
until the second outburst. During the second outburst intensive observation of the source 
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commenced. The greater sensitivity of USA is better able to cover the final decline of the 
source than RXTE. 

In order to proceed with the fitting, we need to select a value for the parameter that 
defines the radial dependence of the viscosity in the disk. We have no strong physical reason 
to prefer one value of over another, so we have chosen ip — 2.8 primarily because it seems 
to work relatively well. We found that close to 3 worked better than of 1 or 2, so we 
chose a ip slightly less than 3 so that equation (62) would still be valid. In the future, as this 
method is applied to other sources it may become clear that there is a preferred value for ip 
which is valid for a class of sources. 

Figure 3 shows the fit of the Green's function (eq. [62], with ip = 2.8) to the first 
outburst. The fit has only two parameters: the viscous timescale to, and the initial time, 
^start- A normalization amplitude is also fitted in Figure 3, but subsequently discarded when 
the function is normalized to unity. Hence the sole nontrivial parameter (beyond our choice 
01-0) in the fit is to- The integral of this measured function is found to approximate the 
shape of the initial rise of the second outburst. The integral has the shape expected from a 
step function turn on; evidently the second outburst begins with something close to a step 
function (see Figure 4). This result is a consistency check of the convolution model for the 
second outburst. In fact, for any input function A{t) whose support is much larger than that 
of the Green's function K(t), the shape of the convolution (eq. [63]) begins with the integral 
of K(t) multiplied by the initial value of A{t). This property of the convolution is a useful 
tool to distinguish the convolution model of the outburst from other possible models. 

Equations (66-67) are used to extract the form of the input function A{t). Figure 5 
verifies that this is the solution by re-folding the solution with the Green's function to derive 
the X-ray luminosity. The solution is smoothed because the convolution with the Green's 
function acts as a low-pass filter. Figure 6 displays the inferred mass Mdit) and its reciprocal 
M^^{t) on a plot of the QPO evolution. The interesting result is that the total disk mass 
is a nearly monotonically decreasing function during the interval when the millihertz QPO 
frequency is increasing. Despite the fact that the computed evolution of the disk mass looks 
very similar to the evolution of the X-ray luminosity, the two curves are not identical. The 
disk mass evolution is obtained from the light curve by the transformation which is effectively 
a shift transformation. 
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7. Discussion and Interpretation of Observable Quantities 

7.1. Evolution of the disk mass and ~ 0.1 Hz QPO frequency 

While Figure 6 is not yet a complete model of the QPO evolution, we argue that it 
represents a very substantial advance over what was previously available. To date there has 
been no model which successfully explains the evolution of this low-frequency QPO. If the 

QPO is some large scale phenomenon in the disk then it should correlate with large scale disk 
properties. It fails completely to correlate with X-ray luminosity in the sense that Lj. rises 
and falls while QPO rises monotonically. The function M'^^(t), while not perfect, moves 
the turnover earlier than that seen for L.^ and hence closer to that observed in the QPO 
frequency history. We believe that this result, combined with a physical model for the QPO 
evolution represents a significant improvement in our understanding of the observed QPO. 
Note that our model included several simplifying assumptions which when corrected may 
account for the remaining differences. 

We may now reconsider the global oscillation model (TOOO,) in which a monotonic 
dependence of MJ"^ on the QPO frequency vq is derived, giving oc z/q (see eqs. [4, 7] 
there). The observed dependence is fiatter than the predicted one; that discrepancy could 
be a result of our assumption of constant conversion efficiency £eff of the mass accretion rate 
(in fact, the surface density) at the inner disk edge into X-ray luminosity. It is possible that 
the inner radius of the disk changes along with the outer disk radius, causing the variation 

of £eflf- 

The global oscillation is a vertical displacement of a large portion of the disk. The exact 
mechanism whereby this modulates X-ray (as well as optical and EUV) emission to produce 
a QPO remains unclear. One possibility is that the modulation is geometrical in nature, 
e.g from varying obscuration of the hot inner portions of the disk. This would require 
the system to be at a fairly high inclination, as found by Wagner et al. (2001) for this 
source. In addition, for this idea be viable the corona would have to be somewhat limited in 
extent. Thus if the relevant disk oscillations were at 10^-Rs (where i?s is the Schwarzschild 
radius) and the corona were confined to less than 10^~^i?s the geometrical model could work. 
Other scenarios could include having the disk brightness respond in some way to the global 
oscillation. 
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7.2. Disk kinematic viscosity estimate and turbulent disk scale 

The best-fit parameter to = 48 d in tlie Green's function from equation (62) (see Fig. 3) 
allows us to derive the kinematic viscosity value because to = 16-Ro/3'^(-Ro) from equation 
(49). Thus the kinematic viscosity is z/ ~ 2 x 10^^ cm^ s^^ since i?o ~ 4 x 10^'' cm (from 
eq. [10] in TOOO for Rout and using m = M/Mq = 7 and the orbital period 4.1 hours). This 
estimate leads to = 6 x 10^ cm if we assume that Vt is of order of a thermal velocity at the 
outer edge of the disk, i.e. Vt ~ 10^ cm s"^. 

8. Discussion and Conclusions 

Here we have presented a detailed mathematical analysis of the diffusion propagation 
model for X-ray outbursts. We investigated the intrinsic properties of the disk density 
evolution equation (5) in a general case. We have analyzed the diffusion models determined 
by the disk kinematic viscosity dependence on the radius. We have demonstrated that the 
Green's (response) function is characterized by the fast rising and exponential decay parts 
(FRED shape) and that it is an intrinsic property for a wide class of diffusion models. 
Similar results were obtained by Sunyaev & Titarchuk (1980) for photon diffusion models. 
Particularly, the light curve of the instant turn-on of a photon source placed in the center of 
a spherical cloud (or the center disk plane) has a typical FRED type shape (Figs. 1-2 there). 

The success in fitting the diffusion model to the observed light curve is very promising. 
We have examined the case where the viscosity is a power law function of position in the 
disk. We note that a dependence of viscosity on surface density is also possible, but this 
will make the problem nonlinear (see Lyubarskii & Shakura 1987 and Lyubarskii 1997 for 
details) . 

We showed in §5 that the steepest descent method can be employed to solve the inverse 
problem (a deconvolution) and extract the input mass accretion function A{t) using the 
observed X-ray luminosity evolution. Such a deconvolution is possible because the inner 
mass accretion rate and the X-ray rate production are related to each other and the timescale 
for viscous diffusion is short enough compared to the timescale of the outburst. Using the 
derived A{t) as the source term and as a measure of the mass lost from the disk into 
the black hole, we derived the mass of the disk as a function of time. This appears to 
inversely correlate with the observed QPO frequency as would be expected from the global 
disk oscillation model, but the correlation is not perfect. 

There are several factors that have been ignored for simplicity in this initial treatment 
which could modify the derived disk mass as a function of time: (1) a changing efficiency of 
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the conversion of mass to at the inner edge of the disk, (2) an effective viscous time (to) 
which varies with time, and (3) an additional sink for disk matter such as a jet. We reahze 
that the model is imperfect but we think that the initial application is encouraging and that 
further study and application to other sources is warranted. 

We conclude by summarizing the main results of the diffusion model. (1) The asymptotic 
form (eq. [62]) is characterized by the exponential rise when the escape time is less than the 
diffusion timescale from the location of mass injection. The same asymptotic form would be 
seen for a semi-infinite medium and thus is determined by the escape boundary condition (in 
our case at i? = /?in) and the location of the mass source. This asymptotic form is insensitive 
to the size of the medium. The second asymptotic form is characterized by an exponential 
decay with the viscous timescale to which is a signature of the bounded medium (in our case, 
a disk) rather than any specific model of the diffusion. (2) The "FRED" shape of the first 
outburst of XTE J1118-I-480 is in good agreement with the diffusion Green's function form 
(eq. [62]), for which only one free parameter to is important. The variation of the viscosity 
index -0 leads to a change of to but the "FRED" shape remains. 

This particular result is very important because it is closely related to the intrinsic 
property of the diffusion Green's function. (3) The value of the parameter to has a physically 
plausible value for the viscous timescale of the disk and ultimately can be applied to the 
calculation of the disk turbulent scale using certain assumptions regarding the disk size and 
the turbulent velocity (4). The rise of the second outburst follows the form expected for 
a step function input, the integral of the Green's function. It is worthwhile to emphasize 
that the shape of this rise is precisely the integral of the Green's function derived from the 
first outburst. (5) The derived mass varies approximately inversely with the QPO frequency. 
(6) Similar X-ray luminosity and QPO evolution would be expected in other sources with 
similar behavior, and in fact is seen in X1354— 644 and GROJ1719— 24 (Brocksopp et al. 
2000). In our model, the rising QPO frequency is a signature of the accretion disk mass 
decreasing as it disappears into the black hole. If QPO frequency varies inversely with disk 
mass throughout the outburst there is a possibility that early detection - before the mass 
reaches its maximum - might find the QPO when its frequency is still decreasing. Both 
the relatively low brightness in the earliest stages of outbursts and the difficulties of prompt 
response would appear to have selected against detection of this effect to date. 
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Fig. 1. — Complete outburst light curve of XTE J1118+480. The time axis is labeled as 
days since 1999 December 23. Diamonds are data from the RXTE/ASM and crosses are 
from USA. The greater sensitivity of USA is reflected in reduced scatter after day 110 where 
the USA data begin. The data have been interpolated and re-gridded so as to give one data 
point per day. 
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Fig. 2. — Comparison of the Green's Function K{t) for various values of ip. The sohd hne is 
the series formulation (eq. [51]) and the dotted line is the combination of asymptotic forms 
(eq. [62]). There is no dotted line in the bottom right panel because equation (62) breaks 
down for ip >3. 
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Fig. 3. — Fit of Green's function (eq. [62]) to initial outburst. For the fit, we have chosen 
^ = 2.8. Other than the trivial amplitude and start time, there is only one free parameter, 
to. 
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Fig. 4. — Plot of complete second outburst showing integral of the best fit Green's function 
from the first outburst as a solid hne. This is the shape that would be expected for a step 
function turn on of the mass transfer causing the second outburst. 
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Fig. 5. — Convolution of derived A{t) with Green's function to get estimated (smoothed) 
Lxit). The solution process integrates and smooths the input like a low pass filter. 
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Fig. 6. — Comparison of the QPO frequency evolution with the derived mass of the disk. The 
points are the QPO frequency according to the left axis (crosses are X-ray determined QPO 
frequencies from Wood et al. (2000) and diamonds are optically determined QPO frequencies 
from Haswell et al. (2000) and J. Patterson & D. Skillman, private commimication). The 
dotted line and the solid line arc Md and M^^ respectively with arbitrary scale factors to 
put them on the same scale as the QPO frequency. The ripple in M^, and Mj^ derives from 
the scatter in the daily X-ray flux measurements, propagated through the inversion process. 
The dashed line is the smoothed included to highlight the shift between and Ma- 



